Regresión: Modelos Estadísticos

Conjunto de Datos: Cheddar Faraway

Abstract

Hemos analizado con las herramientas proporcionadas en el curso de Modelos Estadísticos el conjunto de datos, Cheddar, distribuido en la librería Faraway de R. Para ello hemos utilizado diversas técnicas de regresión lineal.

Introducción

En un estudio de queso Cheddar realizado en el Valle de Latrobe (Victoria, Australia), se estudiaron muestras de queso en las que se analizó su composición química y fueron dadas a probar a distintos sujetos para que valoraran su sabor. Los valores asignados a cada queso son el resultado de combinar las distintas valoraciones.

El DataFrame cheddar de la librería faraway consiste de 30 muestras de queso Cheddar en las que se ha medido el sabor (taste) y las concentraciones de ácido acético (Acetic), ácido sulfhídrico (H2S) y lactosa (Lactic).

Tenenemos un conjunto de datos en el que se recogen observaciones de una cata de quesos, nuestras variables son:

· Taste: una valoración subjetiva de los jueces.

· Acetic: la concentración de ácido acético.

· H2S: la concentración de sulfito de hidrógeno.

· Lactic: concentración de ácido láctico.

Cargamos los datos y enseñamos las primeras observaciones.

Si en nuestro dataset tuviesemos entradas vacías (NA), tenemos varias posibilidades para lidiar con este problema:

• No utilizar/Eliminar las observación que contienen valores.

• No utilizar/Eliminar las variables que contienen las entradas vacías.

• Intentar completar los valores. Existen métodos menos y más sofisitcados:

  • Remplezar con la media, media o moda.

  • Crear una nueva categoría para valores vacíos.

  • Utilizar algún modelo de regresión.

  • Usar un modelo de K-Nearest Neighbors (KNN).

Pero en nuestro caso evaluamos y no hay valores missing.

En nuestro caso, todas las variables son numéricas (cuantitativas). No hay ninguna variable categórica ni cualitativa. Si hubiese, tendríamos que transformarlas en variables binarias. Para ello, deberíamos hacer encoding a variables binarias, el lenguaje de programación R nos permite utilizar as.numeric.

Con estas variables vamos a intentar explicar cómo los valores observados de una variable Y (taste) dependen de los valores de otras variables (Acetic, H2S, Latcic), a través de una relación funcional lineal del tipo Y = f(X). También vamos a intentar predecir el valor de la variable Y para valores no observados de las variables X.

Tenemos 30 observaciones en nuestro dataset. Ahora procedemos a divirlo en el conjunto de train y test. El primero lo utilizaremos para entrenar nuestros modelos y el segundo lo usaremos para cuantificar el error de los modelos.

Para asegurar que sea reproducible utilizamos una semilla, que permite fijar los valores pseudoaleatorios obtenidos en muchas de las funciones utilizadas.

Hacemos un estudio preliminar de nuestras variables. Mostramos un scatter plot de cada variable contrastada con el resto. Esto permite ver a ojo si algún par de variables tiene correlación.

Ahora, utilizamos la función summary de R, la cual nos permite estimar algunos de las características de la distribución del dataset. La siguiente tabla nos muestra los estadísticos más comunes: el mínimo, máximo, mediana, media y el 1er y 3er cuartil.

     taste           Acetic           H2S             Lactic     
 Min.   : 0.70   Min.   :4.477   Min.   : 2.996   Min.   :0.860  
 1st Qu.:13.55   1st Qu.:5.237   1st Qu.: 3.978   1st Qu.:1.250  
 Median :20.95   Median :5.425   Median : 5.329   Median :1.450  
 Mean   :24.53   Mean   :5.498   Mean   : 5.942   Mean   :1.442  
 3rd Qu.:36.70   3rd Qu.:5.883   3rd Qu.: 7.575   3rd Qu.:1.667  
 Max.   :57.20   Max.   :6.458   Max.   :10.199   Max.   :2.010  

Hacemos las gráficas de dispersión entre la variable respuesta taste y las variables predictoras Acetic, H2S, Lactic.

Podemos observar que la que aperentemente guarda una menor relación lineal con taste es la variable Acetic, este hecho será reafirmado más adelante.

Estudio y evaluación del modelo completo.

Intentaremos predecir la variable taste usando el resto de variables. Para empezar, definimos el modelo completo, el cual se usan todas las variables para nuestro modelo lineal múltiple.

\[ Y_i =\beta_0 + \beta_1 x_{i1} + \dots + \beta_{p-1}x_{i(p-1)} + \epsilon_i,\quad i=1,\dots,n \] donde \(Y_i\) es el valor de la variable respuesta para el individuo \(i\)-ésimo, \(\beta_0\) y los \(\beta_j\) son los parámetros \(j = 1, ..., p - 1\), \(x_{ij}\) son los elementos de la matriz de las variables explicativas \(\epsilon_i\) es el término del error aleatorio que suponemos que se distribuye como una \(\mathcal{N}(0,\sigma^2)\), donde \(\sigma^2\) es la varianza que suele ser desconocida.

Resolución mediante matrices

Utilizamos el método de mínimos cuadrados que estima los valores \(\hat{\beta}\) intentando minimizar los errores \(\epsilon\). Como hemos visto en clase, la fórmula que se deduce es:

\[ \hat{\beta} = (X^tX)^{-1}X^tY \] donde \(X\) es una columna de 1’s concatenada con las variables que usamos para predecir.

Por tanto, aproximamos las \(\beta\) modelo lineal completo con los valores de \(\hat{\beta}_0,\dots,\hat{\beta}_3\) con los siguientes valores:

\[ \hat{\beta}_0=-28.8767696,\quad\hat{\beta}_1=0.3277413,\quad\hat{\beta}_2=3.911841,\quad\hat{\beta}_3=19.6705434 \]

Resolución usando librerias de R

Podemos utilizar la función \(lm\), ya programada en R. Definimos el modelo completo, y comprobamos las betas:

\[ taste \sim Acetic \, + \, H2S \, + \, Lactic, \, \, data = cheddar \]

Estudiemos preliminarmente si es un modelo lineal adecuado, para ello comprobaremos las hipótesis estándar del modelo lineal de regresión usando el test de normalidad Shapiro-Wilk. La función \(shapiro.test\) le pasamos por parámetro el residuo/error de cada una de las muestras y nos devuelve un \(p\)-valor.

Observamos que estamos en la hipótesis de que el error nuestro modelo se distribuye de manera normal, ya que el p-valor es \(0.8865 > 0.05\).

Otro aspecto que debemos saber es si la media de los errores será nula o no, para ello usamos una función implementada en R que nos afirma que en nuestro caso esto es así. Respecto a los errores tenemos que saber que se distribuyen con varianza constante, que de manera análoga comprobamos, teniendo p-valores lo suficientemente bajos como para rechazar varianza no constante. De igual manera tenemos que no hay autocorrelación en el modelo completo.

Una vez hemos concluido que aunque estamos en las hipótesis de regresión lineal el modelo completo a pesar de ser el más complejo probablemente da resultados similares a otro más simple.

Call:
lm(formula = taste ~ ., data = cheddar)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.390  -6.612  -1.009   4.908  25.449 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -28.8768    19.7354  -1.463  0.15540   
Acetic        0.3277     4.4598   0.073  0.94198   
H2S           3.9118     1.2484   3.133  0.00425 **
Lactic       19.6705     8.6291   2.280  0.03108 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

s: 10.13 on 26 degrees of freedom
Multiple R-squared: 0.6518,
Adjusted R-squared: 0.6116 
F-statistic: 16.22 on 3 and 26 DF,  p-value: 3.81e-06 

Tras realizar el summary del modelo, una observación muy importante que sacamos es que el p-valor de Acetic es muy grande (0.94) por lo que no va a ser una variable significativa y por tanto con casi toda seguridad no tiene impacto real en el modelo.

Por el contrario el resto de variables sí son significativas por tener p-valor inferior a 0.5, destacando Lactic por tener un beta mayor que H2S, lo que ppuede hacernos llegar a pensar que cobre más importancia en un modelo final. #### Correlaciones y tabla de resultados con el estudio de sus p-valores Usamos el paquete \(GGplot\) de R, el cual nos permite visualizar la correlación y dispersión entre las distintas variables.

Pdemos apreciar que la variable con menor coefeciente de correlación de Pearson es Acetic,mientra que el resto de variables tienen un valor bastante mayor. Cabe destacar que todos tienen un signo positivo, lo que nos dice que en general al aumentar una aumenta la otra.

¿Tiene outliers nuestra muestra?

Para comprobarlo basta realizar el test de Bonferroni sobre nuestro modelo completo:

Concluimos con un nivel \(\alpha\) = 0.05 que no tenemos ningún outlier en nuestra. Lo más cercano a un outlier que tenemos es la observación número 15, que tiene un valor Bonferroni p de 0.17453 (no se acerca a 0.05). Por tanto, no tenemos razones por las que eliminar alguna observación inusual de nuestro conjunto de datos.

Esto se puede comprobar graficamente a través del siguiente gráfico, el cual mide la influencia de cada observación sobre cada una de las betas de nuestro modelo.

Vemos que la que más influye es la antes mencionada observación 15 y por tanto deberíamos considerar si es outlier. Hemos comprobado que sucedería al eliminar esa observación de la muestra, y si bien mejoraba algún p-valor de las hipótesis de regresión, en otros los empeoraba, así que hemos decidido no eliminarla y esperarnos a ver el caso en las seeds que estudiaremos con más detalles si sale de la muestra como observación influyente.

¿Cuál es el mejor modelo?

Como dice el Principio de la Navaja de Ockham, a menudo la explicación más simple es la correcta. Queremos seleccionar predictores que explican los datos de la manera más simple posible, sin disminuir la calidad de las predicciones mucho.

Separacion del dataset en conjuntos de entrenamiento y test (70-30%)

Hemos escogido distintas semillas para estar en condiciones de realizar un estudio más amplio, en la elección de las mismas se ha intentado evitar aquellas que generaban muestras demasiado similares. Las semillas usadas son 1, 1100 y 5 posteriormente se introducirán dos más para afinar en el cálculo de errores.

Consideremos los conjuntos de entrenamiento resultantes de las semillas: train.1(semilla 1), train.2(semilla 1100), train.3(semilla 5).Veaamos que modelos debemos considerar en base a nuestros train.

Método Backward

Partimos del modelo completo estudiado en la sección anterior, eso sí, evaluado en nuestros respectivos train y aplicamos con \(\alpha\) = 0.05, el método de Backward, que consiste en eliminar la variable que menos influya a la predicción. Primero realizamos una iteración explícita del método sobre la primera semilla, y posteriormente se construyen a través de la libreria mixlm de R.

Eliminamos Acetic del modelo debido que su p-valor es >0.05.

Repetimos el proceso con la variable H2S, ya que tiene un p-valor mayor que \(\alpha = 0.05\)

El p-valor es menor que \(\alpha = 0.05\), por lo que hemos concluido, ya que no tenemos suficiente certeza para poder eliminar otra variable.Por tanto, tenemos como resultado que la variable que mejor explica el taste es Lactic. Modelo resultante: taste ~ Lactic, data = cheddar[train.1,].

Por otro lado los modelos backward resultantes por mixlm son:

taste ~ H2S + Lactic, data = cheddar[train.2,] y por otro lado taste ~ H2S + Lactic, data = cheddar[train.3,].

Método Forward

El método Forward consiste en empezar con un modelo de una variable y vamos añadiendo las que más influyan, desarrollaremos el primer modelo de forma explícita y el resto los generaremos con mixlm. De esta manera tenemos:

Actualizamos añadiendo Lactic por tener el menor p-valor.

Con nivel de significación \(\alpha\) = 0.05 este sería nuestro modelo final. Modelo resultante = taste ~ Lactic, data = cheddar[train.1,]
Por otro lado los modelos forward resultantes por mixlm son:

taste ~ H2S + Lactic, data = cheddar[train.2,] y por otro lado taste ~ H2S + Lactic, data = cheddar[train.3,].

Construcción por criterios

En esta subsección trataremos de encontrar un candidato a mejor modelo, construyendo nuestros modelos usando distintos enfoques. Tras aplicar los siguientes criterios a la hora del desarrollo de modelos: \(R^2\) ajustado, Cp de Mallows, Criterio de Informacion de Bayes (BIC), Criterio de Informacion de Akaike (AIC) eligiendo el que nos de un menor o mayor coeficiente según el método (los desarrolos se pueden encontrar en el script), llegamos a las siguientes conclusiones:
solo aparece un modelo nuevo usando el criterio del estadístico \(R^2\) para la primera semilla, taste ~ H2S + Lactic, data = cheddar[train.1,].
Notamos que la combinación de H2S + L aparece en todos nuestros conjuntos de entrenamiento en algún momento, es candidata a ser nuestra mejor elección.

Comparamos los modelos obtenidos hasta ahora en su respectiva muestra de entrenamiento con el modelo completo en ese conjunto de entrenamiento.

Analysis of Variance Table

Model 1: taste ~ Lactic
Model 2: taste ~ Acetic + H2S + Lactic
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     19 1853.5                              
2     17 1410.7  2    442.84 2.6683 0.09821 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Con estos p-valores podemos decir que con un nivel de significación \(\alpha\) ningún modelo es notablemente diferente de su contraparte salvo en el caso de los modelos resultantes en train2 esto puede ser por la cantidad de observaciones influyentes presentes en la muestra, lo trataremos en la siguiente sección.

Diagnostico: Comprobaciones de hipotesis, outliers y observaciones influyentes

En esta sección estudiaremos si nuestros modelos cumplen las condiciones necesarias de un modelo de regresión lineal.

Nuestro enfoque consistirá en un análisis gráfico, acompañado de tests estadísticos en los casos en los que se aprecie una discrepancia notable.

¿Son nuestros modelos, modelos de regresión lineal?: Comprobación de hipótesis.

En la sección 3 se toma un enfoque naïve a la hora de construir los modelos, ya que no hemos estudiado si hay observaciones influyentes, podríamos tener una muestra que no es la adecuada para el estudio de nuestros datos.

Un modelo de regresión lineal debe satisfacer las siguientes hipótesis con nivel de significación \(\alpha\) adecuado:

1. Los errores \(\epsilon_{i}\) tienen distribución normal.
   
2. Los errores \(\epsilon_{i}\) tienen media cero.
   
3. Los errores \(\epsilon_{i}\) tienen varianza constante.
   
4. Los errores \(\epsilon_{i}\) no están correlacionados.

Podemos observar que se verifican a nivel de significación \(\alpha\) = 0.05 se verifican todas las hipótesis de modelo de regresión lineal. En las siguientes gráficas podemos observar como los residuos de los modelos taste ~ H2S + L de los conjuntos train2 y train3 se comportan mejor que cualquiera de los modelos propuestos en la muestra train1

Posteriormente, realizamos un estudio de la colinealidad. Para ello, observamos que el modelo taste ~ Lactic , data=cheddar[train.1,] solo tiene un predictor luego no hay presente ningún tipo colinealidad.

Nuestros valores son muy buenos, entendiendo por bueno VIF < 10, por lo tanto no tenenemos que preocuparnos de una colinealidad grave entre las variables.

Estudio de outliers

A fin de obtener distintos puntos de vista utilizaremos dos métodos, un valor de Bonferroni en un estadístico \(t_{1-\frac{\alpha}{2n};n-p-1}\) y la función outlierTest que utiliza p-valores de Bonferroni obtenidos a través de t-tests.

En los cuatro modelos y bajo los dos criterios no se obtiene ninguna observación que se pueda aceptar como outlier a nivel de significación \(\alpha\) = 0.05

Estudio de observaciones Influyentes

Al igual que en la anterior sección a fin de poder encontrar el criterio que mejor se ajuste a nuestro modelo, en el que hay que tener en cuenta que es una muestra de tamaño reducido, utilizaremos varios criterios, además de distintas técnicas gráficas. Los criterios usados serán:
1.Criterio 1: valores leverage (hii) mayores que \(\frac{2p}{n}\).
2.Criterio 2: valores |DFFITS| son mayores que \(2 \cdot \sqrt{\frac{2p}{n}}\).

3.Criterio 3: valores |DFBETAS| mayores que \(2 \cdot \sqrt{\frac{2p}{n}}\).

4.Criterio 4: InfluencePlot.

En la siguiente tabla observamos las observaciones influyentes obtenidas por cada criterio.

Tomamos la decisión de elegir los elementos que se repitan en varios ya que si eliminásemos los del criterio 3 el resultado sería demasiado pequeño para ser apto para la regresión lineal. Nuestro criterio a seguir es que estén en al menos dos criterios de los restantes.

Estas observaciones son las presentes en la cuarta columna de la tabla y en general son en las primeras en las que nos fijemos a la hora de ver si las tomamos como influyentes o no.

Bajo nuestro criterio quedarían como posibles influyentes las siguientes observaciones:

Ahora es cuando tenemos que determinar si van a ser influyentes o no, esto lo veremos comparando los modelos con las muestras de entrenamiento sin eliminar influyentes y con las muestras modificadas, y evaluar si mejoran o empeoran los p-valores de las hipótesis de los modelos de regresión.

taste ~ + Lactic, data = cheddar[train.1,] vs taste ~ Lactic, data = cheddar[train.1inf,]

Nuestras hipótesis del modelo lineal se ven notablemente mejoradas, tenemos una distribución más parecida a una normal, mejor distribuida en torno a una lineal.
Realizamos un pequeño intercambio en el que nuestra varianza parece menos constante pero se mantiene por encima de todos los niveles de significación habituales.
Tomamos la decisión de actualizar cheddar[train.1,] con cheddar[traininf.1,], que es la versión sin observaciones influyentes. ### taste ~ H2S + Lactic, data = cheddar[train.1,] vs taste ~ H2S + Lactic, data = cheddar[train.1critinf,]

Este es un caso un tanto particular, tenemos unos datos perfectamente distribuidos, pero ya lo estaban antes practicamente. A cambio de eso perdemos un cierto grado de linealidad que hay que tener en cuenta, en este caso no hay una diferencia tan notable con el modelo como en que solo incorporaba a Lactic.

Tomamos la decisión de actualizar cheddar[train.1crit,] con cheddar[train.1critinf,], que es la versión sin observaciones influyentes.

taste ~ H2S + Lactic, data = cheddar[train.2,] vs taste ~ H2S + Lactic, data = cheddar[train.2inf,]

En cierta manera es notable que al eliminar nuestras influyentes y en un dataset tan pequeño no genere a penas diferencia, en todo caso nos mejora todas las hipótesis del modelo lineal, se puede afirmar que es un modelo mejor para trabajar sobre el a priori.

Tomamos la decisión de actualizar cheddar[train.2,] con cheddar[traininf.2,], que es la versión sin observaciones influyentes.

taste ~ H2S + Lactic, data = cheddar[train.3,]

En este modelo, bajo el criterio de elección de observaciones influyentes que elegimos para tratar nuestros datos, este modelo no presentaba ninguna observación influyente, por lo que se mantiene como está.

Errores de Test. Comparacion de Modelos

Nuestro razonamiento para enfrentarnos a esta sección es el siguiente, nos han salido dos modelos posibles y para comprender cual se ajusta mejor a nuestros datos vamos a escoger cinco seeds y evaluar cada modelo en todas ellas, a fin de hacer una media de los errores.

Para esto usamos las tres seeds que hemos utilizado a lo largo del documento, y le añadimos otras dos elegidas al azar. Nótese que evaluar en la primera seed ya está hecho, pues en esta nos salían los dos modelos a considerar. Además en las seeds dos y tres ya está hecho para el modelo con H2S + Lactic, pero hay que repetir el proceso para el otro modelo. De esta manera realizamos los mismos cálculos que los realizados en la parte de diagnóstico, sobre las combinaciones. Esto es, nos aseguramos que cada modelo con todos los train cumpla las hipótesis de normalidad, media de errores nula, homocedasticidad, linealidad y autocorrelación. Después de esto nos planteamos si tiene datos influyentes según las funciones adecuadas, en caso de tenerlos realizamos una prueba. los quitamos del train y vemos si se mejora el p-valor de alguna de las características anteriores y en base a eso decidimos si eliminamos las observaciones influyentes o no.

Con todo eso realizado llegamos a la siguiente tabla, que nos permite asumir como validos todos los casos y calcular sus errores.

Calculamos los errores de cada método como la media de los errores del método aplicado en cada seed. De esta manera, el error medio obtenido es:

Como el error del modelo H2S + Lactic es menor, ese es el modelo que llamaremos final y que se ajustará lo mejor posible a los datos de nuestro fichero.

Conclusión: presentación del modelo final

Finalmente por comparación de los errores sabemos que el modelo final que presentamos es: taste ~ H2S + Lactic , que tiene sentido ya que Acetic no era significativa. otro hecho que lo reafirma es que en la construcción de modelos salieron en todos los casos este modelo, salvo en el método backward y en el método forward de la primera semilla.

Presentamos el summary de nuestro modelo, del que sacaremos más informción en un futuro.

A su vez enseñamos los graficos de este modelo:

Ahora vamos a asegurarnos que verifica las hipótesis para una regresión lineal y además de revisar con outlierTest que no tenemos outliers. También planteamos una hipótesis sobre la nulidad de la media de los errores y vemos que su varianza es constante además de calcularla.

Dicho esto procedemos a la presentación del modelo con sus betas asociados, que si bien se pueden recoger del summary, también los calculamos de forma matricial igual que se hizo con el modelo completo en su momento.

El resultado sería taste ~ -27.591815 +3.9946267 H2S +19.887204 Lactic Nótese que sigue siendo acorde al modelo completo donde el \(\beta\) de Lactic es muy superior en comparación al de H2S.

Presentamos ahora \(R^{2}\) y \(R^{2}_{adj}\), estas como antes se pueden recoger directamente del summary, también las calculamos a partir de los errores y la tabla anova, dando el resultado de 0.6517024 y 0.6259025 respectivamente. Además presentamos el vector de p-valores.

Observamos que todos los p-valores están por debajo de nuestra alpha de referencia.

También presentamos intervalos para las betas de Bonferroni y Scheffé:

Finalmente buscamos una representación de la regresión que tenemos en 3 dimensiones ya que el modelo final consta de una variable respuesta y dos predictoras y vemos el plano de regresión, marcando en rojo las observaciones que peor se ajustan.